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Abstract 

We consider extinction times for a class of birth-death processes commonly found 
in applications, where there is a control parameter which determines whether the pop- 
ulation quickly becomes extinct, or rather persists for a long time. We give an exact 
expression for the discrete case and its asymptotic expansion for large values of the 
population. We have results below the threshold, at the threshold, and above the 
threshold (where there is a quasi-stationary state and the extinction time is very long.) 
We show that the Fokker-Planck approximation is valid only quite near the thresh- 
old. We compare our analytical results to numerical simulations for the SIS epidemic 
model, which is in the class that we treat. This is an interesting example of the delicate 
relationship between discrete and continuum treatments of the same problem. 
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1 Introduction 



Birth-death processes are widely used as a description of processes in physics, chemistry 
[TJI2], population biology and many other areas. These are Markov processes defined on 
states which we label by n = 0, 1, ...,R where R denotes the largest value allowed (which 
could be oo). They are defined by the birth and death rates: 

A n n —> n + 1, 

/i n n— (1) 

The processes we will consider have an absorbing state ('extinction') which we put at 
n = 0; That is, Ao = 0. In what follows we will be mainly concerned with the mean time 
to extinction, r^, i.e., the mean first passage time to the state n = starting at n — k. 
For any Markov process with an absorbing state, extinction will occur as t — > oo with unit 
probability. For example, if n denotes the size of a population of organisms, we seek the 
mean time to biological extinction. 

In this paper we will give exact expressions for the extinction time for a class of birth- 
death processes, and asymptotic expressions for cases where the typical n is large, i.e., in a 
'continuum' limit. We will investigate the validity of a popular approximation, the Fokker- 
Planck or diffusion method £Q . We will see that the Fokker-Planck method gives the correct 
asymptotic continuum behavior of r only in very special circumstances. 

To fix ideas, consider the following two processes taken from the literature of epidemiology 
and population biology. 

• The Susceptible- Infected-Susceptible (SIS) model of epidemiology [I]: Imagine a pop- 
ulation of size N within which n individuals suffer from an infection, and the rest, 
N — n, are susceptible. Suppose the infection rate per contact is A/N, the number of 
contacts is n(N — n) and that the recovery rate is unity (fixing the unit of time). A 
recovered individual immediately becomes susceptible. Then: 

A n = An(l-n/N), 

fi n = n. (2) 

At the deterministic (non-stochastic, continuum) level then there may be a non-zero 
steady state number, n e , of infected individuals, the solution of X n = fi n . In this SIS 
model n e = N(l — 1/A), provided A > 1. This model has a threshold, A = 1, above 
which the infection persists in the contimuum approximation. When the A < 1 the 
infection dies out. In the stochastic model, however, above threshold the number of 
infected individuals remains near n e for a long time (the quasi-stationary state) before 
eventually going extinct [Sj. 
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Figure 1: Regimes for the rates A n ,/i n . a.) Above threshold, b.) Near threshold c.) Below 
threshold. 

• A logistic model from ecology j^j, often called the Verhulst model: This population 
dynamics model assumes a birth rate per individual, B, and unit death rate, per 
individual. In order account for competition for resources, the death rate is assumed 
to increase proportional to n 2 . We write: 

X n = Bn, 

//„ = n+(B-l)n 2 /N, (3) 

defining the carrying capacity, N. At the deterministic level, n e = N provided B > 1. 
In the continuum, for B > 1 the population stabilizes at n e while for B < 1 it goes 
extinct. In the stochastic model there is a quasi-stationary state for B > 1 in which 
the population fluctuates near n e before eventually going extinct. 



These two examples are representative of the class of models that we now consider. In 
both examples there is a large number, N, and we assume that both A and [i involve such a 
number in a special way: 

A n = NX(x), 

/i n = Nj[(x), (4) 
where x = n/N and A, ~p, are smooth functions of x. 

These processes have the following properties: First, we assume that \ n = NX(x) is 
concave downward (or linear) and fi n = N~fi(x) concave upward (or linear). (We will not 
consider the general degenerate case where both functions are linear). Both functions are 
taken to have finite non-zero slopes near n = 0. The processes are most interesting when 
there is a control parameter so that there can be an intersection of the two curves (super- 
threshold), or not (sub-threshold) depending on the parameter. We are also interested in 
the case when the parameter is very near threshold, in a sense that we will define below; see 
Figure 1. 
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We are interested in the mean time to extinction starting at state n. The probabilities, 
ir n (t), for the various states obey the master equation: 

dn n (t)/dt = A n _ivr n „i(t) - (A n + ^ n )iT n {t) + fi n+1 -K n+1 (t) . (5) 

It is an elementary exercise 12] to show that r n obeys: 

— 1 = A n r n+ x — (A n + fi n )r n + /i n r n _i (6) 

We will solve this equation exactly, below, and give expressions for the asymptotic behavior 
of Tfc as N — > oo. In the context of the SIS model there is a very large literature on this 
question IH] • Our general expression agrees with the main results of Ref . |Hj for this 
case. However, in some cases we find different results, see below. 

For N » 1, near the continuum limit, it is tempting to work directly with Eq. (jHJ), and 
try to replace it by a differential equation in x. A naive Taylor expansion gives: 

-1 = f{x)T\x) + ^g 2 {x)T"{x) (7) 

where x = n/N is the continuum variable, T(n/N) = r n , and we define: 

f(x) = X(x) - fi(x) 

9 2 {x) = X(x)+fl{x) (8) 

The smooth functions f(x) and g 2 (x)/2N are, respectively, the drift (sometimes called the 
'force') and the diffusion coefficient. 

The operator on the right-hand side of the first equation above is the adjoint of the 
operator in the Fokker-Planck equation (FPE): PUI2]: 

d t P(x, t) = (l/2N)d xx (g 2 P) - d x {fP) (9) 

This is the evolution equation for the probability transition density of a random walker with 
diffusion coefficient g 2 /2N subject to a drift velocity /. This approach, which is particularly 
popular in the natural sciences, is based on the observation that we can think of the birth- 
death process as a biased random walk in n. From the work of Einstein and Schmoluchowski, 
we can describe a continuum random walk with a diffusion equation, namely 0- The 
extinction time is the first passage time to the origin of the random walker. This approach 
is attractive because the ordinary differential equation in (JSJ) is very easy to solve, and gives 
a tractable formula for T in the limit of large N. Grasman and collaborators used this 
approach for the logistic model [HJ El] , and apparently verified their results by numerical 
calculations. It is also possible to convert Eq. (0) into a stochastic differential equation of 
the Langevin type which has attractive numerical properties. 

However, as we will see below, this series of manipulations only gives the correct asymp- 
totic behavior of r under very special circumstances, namely when the two terms in Eq. (J7|) 
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are comparable in size. In our case, this means that / = A — [i is very small. There are 
several ways in which this can happen. For example, if we are interested in small fluctuations 
near x e = n e /N, then we can almost certainly use a Gaussian approximation and get reliable 
results (though we do not verify this explicitly here). This is what Van Kampen [2] calls the 
'system size expansion', and it is a standard method in applications. 

For the extinction problem we expect that the important n's will be all those between 
n = n e and n = 0. For the force to be small for the whole range we must be near threshold, as 
in Figure (lb). As we will see below, in order for Eq. (J7j) to correctly describe the asymptotic 
behavior we need p n = A n //i n = 1 + (^(l/iV 1 / 3 ^), e > 0, or, equivalently, x e = 0(A r_1 / 3 ~ e ). 
We will discuss the work in Ref. jH] below, and show that the case they treated was, in 
fact, near threshold. This accounts for the the agreement between their FPE results and 
numerical calculations. 

This result is counterintuitive, and demonstrates the delicacy of the continuum limit even 
for the very simple class of processes that we consider here. Explicitly, the continuum FPE is 
useful not for large populations, as one might guess, but only near threshold where the quasi- 
equilibrium population is much smaller than N. The upshot is that a diffusion treatment 
is valid only when the drift is small. In fact, we will go further below: we will show that 
no diffusion description is possible for these processes (except just above threshold) in the 
sense that no diffusion equation can simultaneously give r correctly and also give a good 
approximation to the quasi-steady state. 

In the next section we will state our results. In the next section we give some numerical 
illustrations for the SIS process. In the final section we will discuss the implications of our 
findings for applications. We will relegate the actual computations to the Appendices. 



2 Results 

In this section we summarize our main results. First we briefly discuss the exact solution for 
the mean extinction time as a function of the initial n. Subsequently we present the large 
N asymptotic expansions of the solutions in the superthreshold, threshold and subthreshold 
cases. 



2.1 Extinction time 



1 

T~0 = 0, T R -T R -i = . 



(10) 



The first problem is to solve the second order difference equation © for 1 < n < R — 1 with 
absorption at site n = and a reflecting boundary condition (A# = 0) at site n = R, i.e., 
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This is straightforward: let = r k — t^-i for k — 1 . . . R so that 

-1 = X n a n+ x - p n a n (11) 

for 1 < n < R — 1 with the single boundary condition = Then the solution of this 
first order difference equation for a m , 1 < m < R — 1, is easily written down: 

^ R-m j 

«m = IJPm+i-1 (12) 

A*m . j Pm+j i=1 

where 

Pi = Ai/jUi. (13) 
Then the mean extinction time, r n , is recovered from 



m=l 



By regrouping the product in (|12|) we find the explicit solution with no approximation: 

TO— 1 if, 1 j — 1 



m =l ^ m i=l ^ j=m+l ^ fe=l 



2.2 Expansion for large TV 

Now write A n = NX(n/N) and p n = Np,(n/N) with \(x) and p(:r) uniformly smooth 
functions on [0, r] where r = R/N. We can define p n = p(n/N), where p(x) is a bounded, 
smooth, and non-negative function on [0, r]. We will call the following quantity the 'effective 
potential' 

$(x) = - / logp^K- (16) 
to facilitate comparison to the FPE, below. 
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For large N the products in (|15|) can be estimated by the trapezoid rule of numerical 
analysis. Using the continuous variables z = j/N and y = m/N, 

J] p k = exp J] log p(fc/iV) 

fc=l \fe=l y 

= exp^iV^ logp(w;)du; - ~(logp(0) + logp(z)) + 0(1/ N) 
exp ( JV / log p(w)dw J x (1 + 0(1/ N)) 



J— e-^xft + OOT. (17) 

vp(o)p(*) 



and similarly, 

m— 1 

II - = yfmM ^ x (1 + 0{l/N)) (18) 
i=i pi 

In these estimates and others following, the coefficients of the 0{N~ 13 ) error terms depend 
on the regularity of the rate functions A(x) and fl(x). We proceed under the assumption of 
sufficient smoothness so that the estimates are valid. 

In the following we write x = n/N for the initial point. The large N asymptotic behavior of 
f(x) = r n as given in (|T5)l is very different for the superthreshold, threshold, and subthreshold 

cases: 



Superthreshold case 

When A'(0) > /2'(0) there is a unique 'equilibrium' state n e where A„ e = /i ne , or 
equivalently a unique 'deterministic steady state' x e = n e /N > where A(x e ) = Ji(x e ). 
(In the SIS and logistic models this corresponds, respectively, to the conditions A > 1 
and B > 1.) For the stochastic processes, the extinction time is exponentially large as 
N — > oo. We find that f(x) ~ r„ e = f(x e ) for n = O(N), i.e., x = 0(1). Further: 



2 ^(°M°) x *-"* ( '° X ( 1 + Q(h) (19) 

' N[X(x e )jlf(x e ) - H(x e )fae)] V K N>) 

In this region f(x) is independent of x. For x = 0(1/N) there is a boundary layer 
where f(x) depends on x. However, we have found a simple correction factor which 
allows us to give, for all x G [0,r = R/N), a uniform asymptotic approximation, 

f(x) = (1 - e-^ 10 ^ )*) f(x e ). (20) 

Threshold case 

Here A n < fi n for n > 1, but the derivatives of \{x) and £l(x) at are equal. (In 
the SIS and logistic models this corresponds,respectively, to A = 1 and B = 1.) In 
this critical situation the dominant term in the large N asymptotic expansion of the 
extinction time is oc TV 2 : 



for n/N = x = 0(1). The constant prefactor, C ~ 1.57, is the sum 

(2*Q! 
(2 k k\) 2 (2k + 1 



^ (2k)\ , . 



fe=0 
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Subthreshold case 

Below threshold, A n < fi n for all n > 1, so p(x) < 1 for all x > 0. Moreover, the 
derivative of A (a;) is strictly smaller than the derivative of fl(x) for all a; > 0. (In 
SIS and logistic models this corresponds, respectively, to A < 1 and B < 1.) The 
deterministic version go to extinction rapidly in this case. For the stochastic processes, 
the extinction time is logarithmic in N as N — »• oo: 

f ^= m^m l0tNx+o(1) (23) 

for n/N = x = 0(1). 



3 Numerical estimates and the failure of the Fokker- 
Planck approximation 



In this section we compare our results to numerical calculations, and to the FPE approxima- 
tion for the SIS model. Our numerical results are based on a direct evaluation of the exact 
formula Eq. (|T3|) . 



3.1 Numerical and analytical results for the SIS model: 

As we pointed out above, an example of the family of models which concerns us here is the 
SIS model defined by: 

A„ = An(l- — ), ii n = n. (24) 

X(x) = Ax(l - x), fx(x) = x, (25) 
p(x) = A(l - x) (26) 



For this case: 
so that 
and 



$(x) = - / log A(l - £)d£ = (1 - x) logA(l - x) + x - log A. (27) 
Jo 

It is easy to check that $(x) is a convex function. 

The analytical results for this model based on Eqs. (J19)) . (J21)) . and (J23|) are presented in 



o exact results 
6 - asymptotic approximation 



200 400 600 800 1000 1200 1400 1600 

N 



Figure 2: Comparison of analytical and numerical results for the SIS model in the su- 
perthreshold case. \ogf(x e )/N is plotted as a function of N for A = 1.5,2,3 bottom to 
top. 



the following table. 



A 



r n for n/N = x = 0(1) 



superthreshold A > 1 $'(0) < J^e N ^ K ~ x+1 l^ x (1 + 0(1/ N)) FigHJ 

threshold A = 1 $'(0) = 1.5687^^ + log Nx + 0(1) Fig. H 

subthreshold A < 1 $'(0) > 



^\ogNx + 0(l) 



FigE 



The superthreshold case is compared with numerical simulations in Fig. (j2J). Our formula 
agrees with the results of reference [§| and with the exact results. 

We can go further and test the validity of our error estimate in the first line of the 
table above, and also our treatment of the boundary layer. We define the relative error as 
T/f asy — 1, where f asy is given in Eqs. (JTHJ), (j2UJ). We plot A^ times this quantity in Fig. (jSJ) 
to show that the relative error is of order 1/N and is uniform in x. 

The threshold case is shown in Fig. (j3J); we find good agreement between our asymptotic 
formula and the exact results. For the subthreshold case our formulas do not agree with ||, 
but they do agree with the numerical results. This is shown in Fig. (jSJ). 
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Figure 3: N times the relative error for the superthreshold case of the SIS model as a function 
of x for various N and A = 3. 




Figure 4: Comparison of analytical and numerical results for the SIS model in the threshold 
case. f(x) is plotted as a function of N 1 ^ 2 for various values of x. The symbols are exact 
results of Eq. ()15jl and the dashed lines are guides for the eye. The solid line is the prediction 
for the coefficient of N 1 ^ 2 in Eq. (|21jl . The vertical position of the solid line is arbitrary. 
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Figure 5: Comparison of analytical and numerical results for f(x) in the subthreshold case 
of the SIS model for A = 1/2. f(x) is plotted as a function of log A" for several values of x. 
The symbols are exact results and the dashed lines are guides for the eye. The solid line is 
the prediction of Eq. (|2*3Jl . The vertical position of the solid line is arbitrary. The results 
of [2] are also shown as dotted lines which correspond to x — 1/4,1/2,3/4,1 from top to 
bottom. 



3.2 Fokker- Planck approximation 

The Fokker-Planck approach approximates the finite difference equation JBJ) by the differ- 
ential equation (J7J), solves it, and then extracts the large N asymptotics. This has the 
advantage of producing a tractable and generally useful partial differential equation, Eq. 
Q. However, as we will see, it does not give the correct answers in general for our class of 
processes. In this section we will only be concerned with the superthreshold case. 



3.2.1 Asymptotic estimates 



The solution to Eq. (0) is (neglecting 0{l/N 2 )): 

T [x) = _L [ X e -mV(r)-V(y)) dy +2N 

W Jo 

where 



r e -N(V(z)-V{y)) 
Jy \(z)+p,(z) 



dz dy 



v{x) ~j m+m 

plays the role of the effective potential. 



3 , ^azM dfs - 2 rmu 

Jo 9 2 (0 



(28) 



(29) 
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Figure 6: Comparison of FPE and numerical results for the SIS model in the superthreshold 
case, logf (x e )/N is plotted as a function of N for several values of A. 



The asymptotic behavior, using standard techniques pQ is: 

rW » TM - W x (1 + ° (llN)) ■ (30) 

For the special case of the SIS model we find: 

T(x) ** (X^lp ^-" vw x < x + °^ N » • < 31 ' 

where 

-"<*•> = x * iti + < 32 > 

This expression is quite different from the first line of the table above (which agrees with 
the exact results). In Fig. (jHJ) we show a comparison of Eq. (jHlj) with the exact results. 
It is clear that there is a significant discrepancy for large values of A, that is, far from the 
threshold at A = 1. Note that logr(x e ) /N is plotted in the figure. From the figure we see 
that, for = 1000, T is a factor of about 10 7 smaller than the exact result. 

We have done a similar calculation for the logistic model. The numerical results of Ref. jB] 
are all near to threshold, so that the apparent numerical verification of their FPE calculation 
is only of limited validity. Far above threshold there are similar very large discrepancies 
between T and f . 
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3.2.2 The effective potentials 



We now turn to the source of the problem with the FPE estimate of f. For the superthreshold 
case f in Eq. (j3"Uj) is of the standard form of a relatively slowly varying prefactor multiplied 
by exp(— NV) where V is given in Eq. ()29|) . In the following discussion we will ignore the 
prefactor and consider only the dominant exponential term. The FPE gives a V which is 
not the same as the correct answer, $, Eq. ()16j). 

However, V is close to $ quite near the threshold, in which case force is small over the 
important range of x, namely [0, x e ]. To see this, put X(x) — p(x) = /at(x). Assume that the 
force is, in fact, small: Jn{x) — > appropriately uniformly in x as N — ► oo. Then: 

-*>/*~SN(S)'-KS) ,+ - (33) 

while 

Hence the exponential terms formally agree to a factor of (1+0(1/N)) only if /jv(a^) = X(x) — 
fi(x) < 0(N~ 1 ^ 3 ), that is, the drift must be small over the whole range of x. Equivalently, 
we can write p n = X n / fi n = 1 + 0(l/N l l 3+e ), e > 0, or, x e = O^N- 1 / 3 ^). If this is true then 
the two estimates of f differ by factors of order unity. 



3.2.3 'Corrected' Fokker-Planck equation 



We might be tempted to define a 'corrected' FPE, with an effective potential $ by redefining 
/ and g 2 in Eq. (jHJ). A suitable / and g 2 could be found to do this, but the corrected FPE 
would not be particularly useful. The point of using Eq. Q is to have a unified description of 
the process, equivalent in the continuum limit, to the original master equation. In particular, 
we should be able to describe the quasi-stationary distribution with the same equation. We 
will show that this is not possible. 

A version of the quasi-stationary distribution j3] is obtained by changing the boundary 
condition at the origin to reflecting, that is, we set Ao = 1. Then a stationary distribution 
exists. We can find this by returning to Eq. (jSJ), setting dn n (t) = 0, and solving the equation. 
The result is: 

n^otV/^+i] 

n n = — j (35) 

1 +Er=i rUotv^+i] 

We take the continuum limit by defining p(n/N) = Nir n , and using the method of Eq. 1)17)) . 
We have: 

p{x) x {wW) e ~ NHx){1 + ° {1/N)) (36) 
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We can also solve for the stationary state of Eq. ©, with the requirement that the 
effective potential be $. This gives: 



where C is a constant. We can also recalculate T, as in Eq. (|3()|) . If we compare to Eq. ()19j1 
we see that while we have the correct effective potential (by construction), we do not get the 
correct prefactor, so that T is inconsistent with f. 



4 Summary and Discussion 

We have described two sorts of results in this paper. On the one hand, we have shown how 
to generalize previous work on the SIS model IB1 03 E] to general birth-death processes, 
and to the threshold and subthreshold cases. We think that we have corrected an error in 
Ref. P for the subthreshold case. 

Our treatment of the FPE should be viewed mainly as a warning about the subtlety of 
the relationship between discrete and continuum approaches. In cases where the continuum 
equation, Eq. © 'should' work, it fails if the drift is too large. This is reflected in the fact 
that the exact effective potential, <3> = — f [log(l + f/g 2 ) — log(l — f /g 2 )], is different from 
the FPE expression, V = —2 f[f/g 2 }- They are the same for small /. For practical purposes, 
if the process in not immediately above threshold, it is preferable to use Eq. (|19|) or even 
the exact expression, Eq. (fTKJ) . 

However, for a problem with a wide separation of scales, this option is often not available. 
For example, in work on modeling calcium waves in cells [H], the underlying processes are 
often too complex to allow a practical exact calculation. One method to circumvent this 
difficulty is to use a Langevin equation, replacing some of the rapid processes by noise terms. 
In practice, however, there is some evidence that such a method gives acceptable results only 
near an appropriately defined threshold [T2] . 

We have tried to find a heuristic argument for the failure of the FPE for the extinction 
time. We speculate that in the large N limit, the corresponding stochastic process will have 
discontinuous sample paths since the fluctuations that lead to extinction probably occur over 
times that do not scale with N. In this case, no diffusion approximation is possible. 

Acknowledgements: This work was supported in part by NSF Award DMS-0244419. 



P(x) oc e- N ^ x) /g 2 



(37) 



Comparing these two equations we see that we must set: 




(38) 
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Appendix 



In this appendix we give the details of the large N expansion of Eq. (|15|) which we repeat 
here for convenience: 



771—1 R 1 j — 1 



m=l ^ m i=l " i=m+l ^ fe=l 

In the following sections we discuss the superthreshold, threshold, and subthreshold cases. 
The major computations are of the term A-m in Eq. . 



A Superthreshold 



In this situation is convex with a quadratic minimum at x e , the unique solution of 

\{x e ) = fl(x e ). Thus, invoking standard integral estimates, 

E = E 1 , J_ r/ e-^x(l + Q(l/iV)) (40) 

x ( 1 + (1/N)). 



For the next step we use the fact that if h is a smooth function with h(x e ) ^ 0, then for 
s <C x e , 

^(Oe-^de = h(x e ) ^J^) e~ N ^ x (1 + 0(1/N)) (41) 

Hence, recalling s = m/N, and using n = O(N), $ (0) < 0, $(0) = and the smoothness of 
p to compute the geometric sum after expanding <&(m/N) around 0, 



m=l 



2, v^L x(1+0(1/jV)) ^ (42) 



W(* e ) ^(xe)v^) l-e*'(0) 

The last expression is exponentially large in N so Ylm=i which we estimate in the 
subthreshold case, below, is completely negligible. Eliminating p, then, we find 



27rA(0)/i(0) e -N*( Xc ) n(um \ 

Tn - V N[X(x e )H'(x e )-X'(x e )fl(x e )] A(0)-/i(0) X [ + U[ ,JS)h { S) 
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independent of n for n = O(N). The boundary layer correction follows simply from a finite 
geometric series approximation to the sum in the penultimate line in (J42I) . 



B Threshold 



Here the potential function is convex with vanishing derivative but nonvanishing cur- 
vature at x = 0. Referring to the terms in the exact solution (jSHJ), invoking the trapezoid 
approximation, and putting s = m/N, 



TA m = V f V^MO e -NW»-WO) dz x (1 + 0(1/JV)) 

~y ~,Jm/N\/p(z)ll(z) 



1 X J m/N y/p(z)fJ,(z 



= M f ds f dz e -*(*(*)-*(*)) x (i + o(i /iV )). (44) 

Jo Js \J\(z)n(z) 

Most of the contribution to the integral comes from a neighborhood of the origin, so we 
proceed making standard integral approximations. Expand p(s), X(z) and using 

A(0) = /x(0) = $(0) = $'(0) = 0: 



Y J A m = N^y F p^= ds dz- x (1 + 0(1/JV)). (45) 



Switching the integrals, 



n /^7?V\ /t „-7V<S>"(0).z 2 / 2 /•roinfox) . 



\AW(o) io 

and changing variables 



/ £ dw / e « 2 / 2rfM x (i +0 (i/jv)). (47) 



$"(0) ^'(0)^(0) io 
Then because N is large we may modify the limits of the integrals 



y Am= JL_ VW) r d v < —— [ V due u2 ' 2 x (1 + 0(1/7V)) . (48) 

Integrate e" 2//2 using its Taylor expansion and use the formula for the even moments of a 
gaussian to obtain 



m=l 



$"(0) v/A'(0)/i'(0) ^(2^!)2(2A; + 
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The sum above converges to C ~ 1.5687, so 



m=l 

We show in the next section that 



2$"(0) V A'(0)/x'(0) W 



E— = ^>gn + 0(l). 



Thus, 



C Subthreshold case 



Here the potential function is a convex function with positive derivative at 0. 
= v 7 *^) V P(m/AQ J "" (1+0(1/JV)) 

Ft 

= pim/N)-^ £ w#wf=7ff (l + 0(1/iV)) 

j=m +l Ny/\(j/N)(J,{j/N) 

Adding and subtracting 

p(m/N) 



^ n(m. /mi-fn+l/2 



j=m+ 



1 Ny/\(m/N)ji(m/N) 
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and recalling that in subthreshold case p(x) < 1 uniformly for x G [0,r], we have 

A m = p(m/N)-^ V ( - P{m/N)3 P{m/N)3 ) (55) 



p(m/N) j - m+1 / 2 
1 Ny/\(m/N)p,(m/N) 



j=m+l 



p{m/N) 1 



+OQ./N) 



1 - p{m/N) Npim/N) 

S v ' 

Cm 

= £ m + C m + C»(1/7V). 

Now we evaluate 5^m=i -B m and 5^m=i C m separately as n = Nx — > oo. Changing vari- 
ables in the integral in S m and using the facts that for small s, \/A(s)/i(s) = ^/A(0)/i(0) s + 
0(s 2 ), or y/\(s)ji(s) = y/\(0)fi(0) s x (1 + 0(s)), and p(s) = 0(1), we find 

a = rt*)"* +1/2 _ rt<* +1/2 w m 

Jo W\(Z + s)p(Z + s) V\(sjp{s) ) 



r—s 



r- 



0(1) x / e^iogpW £ m 



< 



0(1) X 



pr/s-l 

= 0(1) x / e {sa)z dz 

Jo 1 + z 

S 0(1) x / ' e {sa)z zdz 
Jo 



[sa] 2 



where we denote a = N log p(s). Then recalling s = m/N, we conclude 

n n 

£ B m ^ 0(1) x £ = 0(1) (57) 

r^l m^l 171 ( lo § P( m / N )) 

as iV — > oo (which means n — > oo as well). 



On the other hand the C m terms may be written 

n n 

111 

m=l 

1 



^ C '=^^m^m (58) 

m=l ' 



where we define 

V(-) = P{m,N) (59) 
V[ N } l-p(m/NY 1 j 

Then subtract and add the 'divergent part' of the sum: 

Vr V r ^ 1 i V ^ (0) (m) 

^ m ^ [ Nh(m/N) mu'(0) J ^mu'(O) 1 J 

m=l m=l i \ i / m=l v ' 

rx i!«)«/i'(o)-i?(o)/2«) , ij(o), , , . 



y ^(o)/i(0 //'(o) 



where 7 = 0.5772.. is Euler's constant. 



Finally, 



N' 



n ^ n j 

= (61) 



O(l) 

log (n) + 0(l). 



/i'(O) 

Putting these calculations together we conclude 
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